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Abstract 

n 

We apply the linear response theory developed in 1] to analyze how a periodic signal of weak 
amplitude, superimposed upon a chaotic background, is transmitted in a network of non linearly 
interacting units. We numerically compute the complex susceptibility and show the existence of 
specific poles (stable resonances) corresponding to the response to perturbations transverse to the 
attractor. Contrary to the poles of correlation functions they depend on the pair emitting/receiving 
units. This dynamic differentiation, induced by non linearities, exhibits the different ability that 
units have to transmit a signal in this network. 

PACS numbers: 89.75.Hc,02.70.-c,05.10.-a,05.45.-m 



I. INTRODUCTION 

Currently, there is considerable research activity in network dynamics. This is clearly 
motivated by the wide expansion of communication media (mobile phones, Internet, mul- 
timedia, etc.), but also by the growing interest in network modeling of biological processes 
(neural networks, genetic networks, ecological networks ...). A large part of these studies 
focuses on topological properties of the underlying graph. However, in many cases, the 
nodes of the networks are units behaving in a non linear way. For example, in a communi- 
cation network a relay regenerates (amplifies) weak signals, but it has a finite capacity and 
saturates if too many signals arrive simultaneously. A neuron has a non linear response to 
an input current, a gene expression is determined by a non linear function of the regula- 
tory proteins concentration, etc.. These non-linearities might modify the network abilities 
in a drastic way. For example, a relay may have a high graph connectivity ("hub"), but 
the dynamics drives it close to its saturation point, so that it has a weak reactivity to the 
changes in the inputs coming from the other units and a poor capacity to transmit informa- 
tion. Consequently, the information is transmitted via other units, possibly weaker links, 
and, in this regime, these units become temporary "hubs" though they may have a low 
graph connectivity, while the main hub is decongested. In biological networks similar effects 
may arise. For example, the capacity of a neuron to transmit a specific excitation strongly 
depends on its state, determined itself by the overall currents coming from afferent neurons. 

This suggests us that the mere study of the graph topological structure of a network 
with non linear units is not sufficient to capture the full dynamical behavior. However, 
there are relatively few studies which analyze the joint effect of topology of the network 
and non-linearity. Nevertheless, these networks are dynamical systems with a large number 
of degrees of freedom, and so dynamical systems theory and statistical mechanics provide a 
powerful framework to state problems in a well-defined way and to propose solutions. 

In this paper, we analyze the following situation. We consider a network composed by a 
set of N units receiving and transmitting signals. At each time step t the unit i receives a 
bench of signals coming from each unit connected to it, and it emits, at time t + 1, a signal 
which is a sigmoid function of the global input [see eq. (J2J)]. In the model studied below, the 
global dynamics has generically a chaotic attractor, provided that the non linearity of the 



transfer function is sufficiently large [see Section |TTj. In spite of the presence of chaos it is 
possible to analyze how a periodic signal of weak amplitude, superimposed upon a chaotic 
background, is transmitted in the network. However, as discussed above, this analysis 
requires the consideration of the network structure as well as nonlinear effects. 

The main tool we use for this investigation is the linear response theory developed by 
D. Ruelle |l| for hyperbolic dynamical systems (e.g. dissipative systems with a chaotic 
attractor) in a non equilibrium steady state. This theory allows us to compute explicitly 
the variation of the average value of a generic observable, induced by a time dependent 
signal of weak amplitude. Indeed, provided that the amplitude of the signal is sufficiently 
small (but finite), this variation is a linear function of the signal and a linear response 
operator is explicitly given in terms of the dynamic evolution. In our case, this operator 
has a simple expression (see eq. ©)• The effects of a periodic signal emitted by a unit 
on a receiving unit is characterized by the Fourier transform of the linear response, called 
susceptibility in the sequel (see section ITVJ). This gives us a frequency response curve (see 
Fig. Q) exhibiting resonances peaks. These resonances corresponds to complex poles for 
the analytic continuation of the susceptibility in the complex plane. They have a nice 
interpretation in Ruelle theory. 

Indeed, in this theory, the linear response operator is the sum of two contributions. There 
is a regular term, corresponding to the response to perturbations "parallel" to the attractor 
(more precisely locally projected along the unstable manifold). This term is actually a 
correlation function j2| and, consequently, it obeys classical relations such as the Fluctuation- 
dissipation theorem. The poles of its Fourier transform are called Ruelle-Pollicott resonances 
[3( or "unstable" poles. They give the rate of mixing of the chaotic system or equivalently, 
the relaxation rate to equilibrium for a perturbation "on" the attractor. These poles are 
independent of the observable. Therefore, in our case, they are independent of the pair 
emitting/receiving unit (see Fig|2j). When focusing on the response to real frequency one 
observes therefore resonance peaks common to all pair of units, and these peaks are also 
present in the Fourier spectrum of the corresponding correlation function. 

The second term corresponds to the response to perturbations locally projected along 
stable manifolds, namely transverse to the attractor. Therefore, it exists only in the 
dissipative case. It does not obey fluctuation-dissipation theorem and has drastically 
different properties than the first term. In particular its poles ( "stable" poles) are expected 



to be distinct from the unstable poles. In this paper, we indeed exhibit such stable poles. 
To the best of our knowledge, this is the first example where these poles are explicitly 
exhibited, though their existence was theoretically proved. Moreover, we show numerically 
that the stable poles depend on the pair emitting/receiving unit (see Fig. |3J). When 
focusing on the response to real frequency one observes therefore specific resonances peaks 
(see Fig. Q). This shows that a unit receiving a periodic signal emitted from another unit 
may respond in a specific way to this signal, the amplitude depending both on the signal 
frequency and on the emitting unit. Note that according to the discussion above this effect 
cannot be observed by studying correlation functions. 

The paper is organized as follows. In section |H] we introduce the model and discuss its 
properties. The section IHII recall briefly the main results of Ruelle linear response theory 
used in this paper. An explicit computation of the linear response is performed. It shows 
the explicit contributions of the network topology and of the non linearity in a signal propa- 
gation. In section HVl we compute numerically the frequency response curve and discuss the 
different resonance peaks. The poles of the complex susceptibility for a few pairs of units 
are computed and compared in the section Our main conclusions are then drawn. 

II. MODEL 



Consider the fo 

Networks [see 



lowing dynamical system, originally proposed in the context of Neural 
and references therein] . The output signal is a function of the weighted 
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sum of the signals arriving at % at time t and is given by: 

N 

Ui (t + i) = Y,JijfMt)) (!) 

The weights J^'s may be positive (excitatory), negative (inhibitory) or zero (no direct link 
between j and i). They are in general non symmetric (Jy 7^ Jji). Thus, the matrix of 
weights, J, defines an oriented graph such that there is a link from j to i if and only if 
Jij 7^ 0. The global dynamics can also be written as: 

u(t+l) = G[u(*)]=Jf(u(t)), (2) 

where u(t) = {ui(t)} i=1 and where we used the notation f(u(t)) = {f(ui(t))} i=1 . Consider 
now the case where the nonlinear transfer function / is a sigmoid, [e.g. f(x) = tanh(gx)], 



where the parameter g controls the non linearity. In terms of input/output ratio, a unit 
amplifies weak signals (if g > 1), but with a limited capacity: / "saturates" if the local field 
is too strong, and the variations of the output signal are all the weaker as the local field is 
big. Thus, the capacity of i to retransmit a signal emitted from k does not only depend on 
the weight J^ but also on the state of saturation of i when it receives the signal coming from 
k. Note also that the Jacobian matrix DG(u) writes DGjj(u) = Jijf'{uj) where /' is the 
derivative of /. Therefore, the volume variation is proportional to Y\f = if'{ui). Therefore, 
in this model, the dynamical contraction is closely related to the saturation of the sigmoid 
transfer function. 

In order to emphasize the effects of the nonlinearity and minimize the effect of the 
network topology, one may assume that the network is fully connected and that the Jy's 
are drawn randomly with respect to some smooth distribution (uniform, Gaussian, etc . . . ). 
As an example, one may fix the average value [Jij\ = and the variance [</?■] = -^ (to ensure 
the correct normalization of the local field with the size N). This example is interesting 
because the system (J2J) exhibits a wide variety of dynamical regimes (static, periodic, quasi 
periodic, chaotic). More precisely, it has been shown in 5} that it generically exhibits a 
transition to chaos by quasi periodicity when g increases. Note that the same transition 
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occurs if the network is sparse j^| with K > 2 neighbors (K can be random) chosen at 
random, provided the variance of the J^-'s scales like j^. However, we do not address this 
case in this paper since we want to minimize the effect of the network structure. Note 



also that this type of transfer function al 



coexist. It has been indeed shown in 



0,3 



ows dynamical regimes where several attractors 



that, adding a threshold 9 to the local field, 
there exists a region in the parameter space g, 8 where two attractors coexist. This region 
can be analytically computed. However, in the present paper, the parameters are located in 
a region where there is only one attractor and all initial conditions converge to this attractor. 

Let us now assume that the non linearity is large enough so that the global dynamics has 
a chaotic attractor (with all Lyapunov exponents bounded away from zero and at least one 
positive Lyapunov exponent). 

We now add a signal of small amplitude £(t) to the output of some units. Then the 



evolution of the perturbed system, denoted by u, is given by : 

u(t + 1) = G [u(t)] + £(t) = G [u(i)] . (3) 

Note that the formalism introduced below accommodates the generalization where £(t) 
depends also on u(t), but we do not consider this case here. 

We want to investigate the capacity of the network to transmit signal £(t) superimposed 
upon the chaotic background. This is a complex problem since after a few time steps the total 
signal arriving at time t at k includes the sum of contributions corresponding to different 
paths followed by £, with different time delays. Moreover, along a path the signal can be 
damped if / saturates (/' < 1), or amplified (/' > 1). Finally, the dynamics being chaotic, 
after a sufficiently long time the signal is distorted by the nonlinearities and scrambled by 
mixing. 

To tackle this problem we analyze how the difference u(t) — u(t) between the perturbed 
and unperturbed dynamics behaves on average as a function of £(£). When £(t) is small 
enough, and in spite of the initial condition sensitivity intrinsic to chaotic systems, it can 
be shown that this difference is a linear functional of £(t). This is the content of the linear 
response theory developed by D. Ruelle |l| for chaotic and dissipative [10] system. Some 
aspects of this theory are briefly recalled in the next Section. 

III. LINEAR RESPONSE THEORY. 

The unperturbed dynamical system (J2J) has a strange (chaotic) attractor for sufficiently 
large g. Usually, strange attractors carry a natural probability measure called the Sinai- 
Ruelle-Bowen (SRB) measure pj. If one prepares the system (J2J) with an initial macrostate 
p having a uniform density (i.e. p(du) = du), corresponding to selecting typical initial con- 
ditions, then, provided that the limit exists, the SRB measure is the asymptotic macrostate 
p = lim t ^ +00 G*/i where G l p is the image of p under the t-th iterate of G. The SRB measure 
has several remarkable features which make it "natural" |ljj. One of its most important 
property for practical purpose is the following: If A is some observable (a smooth function 
of u), its average with respect to p, 

<A>= [ A(u)p(du) (4) 



is equal to the time average along typical trajectories. This means that "ensemble average" 
and time average are equivalently for typical trajectories. This is especially useful for 
numerical computations (see next Section). 

Applying a time dependent perturbation £(£) to the system induces time dependent 
changes in the statistical averages. More precisely, the natural extension of the SRB mea- 
sure defined above is a time dependent SRB measure p t . It is given by the (weak) limit 
lim„ + „ &... G-„. The corresponding average wil. be denoted by < >, 

It has been established in |l| that a linear response theory exists for uniformly hyperbolic 



diffeomorphism T t 12]. In our framework, this means that, provided that £(t) is sufficiently 
small, and for any smooth observable A, the variation < A > t — < A > is proportional to 
£(£) up to small non linear corrections. In other words, pt is differentiable with respect to 
the perturbation. The derivative is called the linear response. 

The theory developed by Ruelle allows one to compute the linear response, for general 
perturbations depending both on time t and state u, and for a general observable A. In our 
context, however, where the considered observables are simply the variables of systems (J2J)- 
(J3J), the linear response has a simple form, which can be written as: 

oo 

<u> t -<u>= £ x ( T )£( t -r-l) (5) 

r= — oo 

where x( T ) represents the averaged Jacobian matrix: 

X(r) =< DG T (u) >, (6) 

for r > 0. Otherwise x( T ) = (which is consistent with the requirement of causality). 

A remarkable consequence of Ruelle theory is that x( T ) is a bounded function for all 
r > 0. In particular, it does not diverge exponentially fast, despite the presence of a positive 
Lyapunov exponent. As discussed below, this is essentially a consequence of exponential 
mixing. 

In what concerns network dynamics, equation (J3J) is interpreted as giving the average 
response of unit % of the system when the network is submitted to weak signal £(t). In 
particular it is seen that if only one unit j is perturbed at time t = — 1 by a kick of 
amplitude e [that is £(t) = eej8(t + 1) with the Kroenecker symbol 5 and the j-th unit 



vector ej], then exijif) gives precisely the mean response of unit i at time t. This suggests 
to define the susceptibility of the network as the Fourier transform of Xij (f), namely: 

oo 

X(u) = £ X (t)e^ (7) 

t=— oo 

This matrix function will be numerically computed and studied in the next Section. We 
conclude the present Section by analyzing further the structure of Xij( T ) i n the case of 
dynamical system (JTJ. Here one can decompose Xijij) as : 

x«(r)= e nv^/n/'K-^-i))), (8) 

The sum holds on each possible paths 7y(r), of length t, connecting the unit A; = j to the 
unit k T = i, in r steps. One remarks that each path is weighted by the product of a topological 
contribution depending only on the weight Jy and a dynamical contribution. Since, in the 
kind of systems we consider, functions / are sigmoids, the weight of a path 7y(r) depends 
crucially on the state of saturation of the units ko, ■ ■ ■ , k T -\ at times 0, . . . , r — 1. Especially, 
if f'iuk^il — 1)) > 1 a signal is amplified while it is damped if / / (wjfe J _ 1 (/ — 1)) < 1- Thus, 
though a signal has many possibilities for going from j to i in r time steps, some paths 
may be "better" than some others, in the sense that their contribution to Xij{ T ) is higher. 
Therefore eq. (jHJ) underlines a key point discussed in the introduction. The analysis of signal 
transmission in a coupled network of dynamical units requires to consider both the topology 
of the interaction graph and the nonlinear dynamical regime of the system. 

IV. COMPLEX SUSCEPTIBILITY. 

One can decompose the response function (jHJ) into two terms. The first one is obtained by 
locally projecting the Jacobian matrix on the unstable directions of the tangent space. This 
term will be named the "unstable" response function. It corresponds to linear response of 
the system to perturbations locally parallel to the local unstable manifold (roughly speaking 
perturbations "on" the attractor). One can show that the linear response associated with 
this type of perturbation is in fact a correlation function, as found in standard fluctuation- 
dissipation theorems [l|. Hence, as usual for correlation functions of a chaotic system, it 
decays exponentially (because of mixing) and the decay rates are associated with the poles 
of its Fourier transform. More precisely, these exponential decay rates correspond to the 
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imaginary part of the complex poles of the unstable part of the susceptibility (JBJ). Thus they 
will be called "unstable" poles. More generally, it can be shown that these poles are also the 
eigenvalues of the operator governing the time-evolution of the probability densities (which 
we denoted above as G 4 //), the so-called Perron- Frobenius operator |3|]. Therefore, these 
poles, whose signatures are visible in the peaks of the modulus of the correlation functions, 
do not depend on the observable, though some residues may accidentally vanish for a given 
observable. 

The second term [l3| is obtained by locally projecting the Jacobian matrix on the stable 
directions of the tangent space. It corresponds to the response to perturbations locally 
parallel to the local stable manifold (namely transverse to the attractor). Therefore, it is 
exponentially damped by the dynamical contraction. [Note that, according to the specific 
form of the Jacobian matrix, this contraction is, in our case, mainly due to the saturation of 
the sigmoid transfer function]. The corresponding exponential decay rates are given by the 
complex poles ( "stable" poles) of the stable part of the complex susceptibility. But here the 
poles depend a priori on the observable. One can easily figures this out if one decomposes 
the stable tangent space of a point in the orthogonal basis of Oseledec modes (directions 
associated to each of the negative Lyapunov exponent). The projection of the i-th canonical 
basis vector on the fe-th Oseledec mode depends on i and k. This dependence persists even 
if one takes an average along the trajectory, as in (JBJ). 

Hence, both stable and unstable terms are exponentially damped, ensuring the con- 
vergence of the series (jHJ, but for completely different reasons. Moreover, the stable and 
unstable part of the linear response have drastically different properties. As a matter of fact, 
the stable part is not a correlation function and it does not obey the fluctuation-dissipation 
theorem. In particular, the unstable poles and stable poles are expected to be distinct. 
In this paper, we give for the first time an evidence of this theoretically predicted effect. 
Moreover, we show that the stable poles indeed allow to distinguish the units in their 
capacity to transmit a signal. 

For this we first numerically compute the susceptibility (j7j) for real values of u. The com- 
putation is based on the following remark. Let us consider perturbations ^'(t) = ee,- cos(ut) 
and £^(t) = —eejsm(uit) and let vS 1 ', u^ 2 ^ denote the variables of the corresponding per- 



turbed systems: 



,( fc ) ( 



;{k) 



u w (£+l) = G u w {t) +e k >{t) (A: = 1,2) (9) 



Then it follows from JSJ) that : 



W ^, - ^ 7,. «0 -i-*/> „( 2 ) ^. - <- «, ^ = ^V > -v..C^ ( o-^(*- r - 1 ) 



;< uj J > t - < Mi >) + i(< uj ' > t - < Ui >) = e^Xij 



r e 



= e&Me-"^ (10) 

Note that the time dependent average response to periodic perturbation is therefore periodic. 
The linear response at time t is an infinite sum corresponding to contributions of time 
delayed signals following different paths. Since the signal is sinusoidal the terms in this sum 
may interfere in a constructive way [but exponential damping prevent the series to diverge, 
ensuring the existence of a linear response] . 

Since Xij(u>) is independent of t, then it is equal (for u ^ 0) to the time average 

*«(") = lim ^ E e*" (w) [< u,^ >t +* < 4 2) >t] (11) 

t^oo Te i=0 

The time-dependent averages < u\ > t involve an average over initial conditions in the 
distant past. One can reasonably assume that the above average over t makes the average 
over the initial conditions unnecessary. Then one may write : 

*«(")= lim ^rj:^^[uf\t)^iut\t)] (12) 

where the u\ (t) (k = 1,2) are obtained by iterating maps (JHJ). This provides a straight- 
forward way to compute the susceptibility, where most of the computing time goes into 
computing the orbits u^ k \t). 

As an example, we performed a numerical computation of the dynamical system (J2J) where 
we take a fixed realization of Jj/s, with N = 8 units. There is a quasi-periodic transition to 
chaos as g increases. The system is studied for g = 3.5 corresponding to a positive Lyapunov 
exponent Ai = 0.158, while the second one is A2 = —0.183. The system is therefore weakly 
hyperbolic (all Lyapunov exponents bounded away from 0). 

The function x(co), the Fourier transform of the matrix (jSJ), has been computed with a 
resolution 5u = ^g ~ 1.53 x 10~ 3 . The average is done with 26214400 samples. We did 
several runs where we varied the length T of the time average in (fT2*|) . We checked that the 
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global structure is the same. In particular the amplitude of the susceptibility |x(w)| does 
not depend on T (see note [12]). Also the fluctuations decrease like -4= according to the 
central limit theorem. 

In Fig. E]we have plotted the modulus of the susceptibilities X33, X45, X71. Comparing 
these functions, one remarks that there are thin peaks essentially located at the same fre- 
quencies, with different heights. Moreover, these frequencies are harmonics of a fundamental 
frequency (u!q ~ 0.166). This is expected from the frequency locking in the quasi-periodic 
transition preceding chaos. Some of these frequencies are also present in the Fourier spectrum 
of the correlation functions but with a smaller amplitude and some peaks are indistinguish- 
able from the background. Instead, all harmonic peaks are revealed in the susceptibility 
spectrum. 

But we also note that for many peaks, the width varies strongly from a pair ij to another. 
This means that the resonance strength depends on which unit is excited and which unit 
responds. In particular, some peaks are very thin, corresponding to an accurate resonance 
while some others are broad. In terms of poles, this means that the imaginary part are 
distinct and consequently the corresponding poles are different [see next section]. Finally 
there are additional peaks strongly dependent on the pair ij. 

Thus, a simple glance to Fig. [T] tells us that the frequency response of a unit i to the 
excitation emitted by a unit j strongly depends on the pair i,j. As discussed above, and 
numerically shown below, this difference comes from the stable part of the linear response. 
Consequently, the specificity of the response is revealed only if one consider perturbations 
transverse to the attractor. [Note that, generically, the signal is a perturbation having local 
projections both on local stable and unstable spaces.] 

V. UNSTABLE AND STABLE POLES. 

Resonances correspond to poles in the complex plane. As a matter of fact, the position 
of the maximum of the peak corresponds to the real part of the pole, its width is related 
to its imaginary part, and the value of the maximum is related to the residue. From this 
observation, we developed an algorithm to estimate the residue width and locations of the 
poles. Let uq = uj r + iuJi be a simple pole of x an d A its residue. If one multiplies x by 
a phase factor e 1 ^ then the real and imaginary part rotate continuously, without changing 
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FIG. 1: (Color on line)Modulus of the susceptibilities X33 (red),X45 (blue) ,£71 (green). 

the modulus. If the pole is close enough to the real axis then there exists a phase ifi such 
that, on the real axis, the real part has a characteristic Lorentzian shape symmetric with 
respect to u r while the imaginary part is antisymmetric. Then a nonlinear curve fitting 
allows us to determine A,u r ,Ui. Once a local analysis has roughly determined the poles, a 
global nonlinear fit (Levenberg-Marquardt j9|) allows us to localize the poles with a better 
accuracy. 

In fig. |21 we have plotted the real and imaginary part of the poles of several correlation 
functions. One notices that all pair of units have poles at the same value of u, within the 
error bars. We have also plotted in Fig. El the modulus of the susceptibilities X33;X36?X63 
(left column) and the corresponding poles (right column) with the poles of the correlation 
functions. As expected from Fig. Q we observe common poles (unstable poles) but also 
distinct poles (stable poles) that, moreover, strongly depend on the pair receiving /emitting 
unit. 

Finally, note that some poles are very close to the real axis. Since their imaginary part 
is related to the coherence time of the response to a kick, this tells us that the response to 
a pulse may subsist on quite a bit long times, though the underlying dynamics is chaotic. 
[Recall however that the linear response measures variations of average value of observables] . 
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FIG. 2: (Color on line) Poles of several correlation functions. 
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FIG. 3: (Color on line) Left column: Susceptibilities X33>X36>X63 and reconstruction by the non- 
linear fitting algorithm (NLF) used to compute the poles. Right column: Poles of susceptibility 
(red squares) and poles of correlations (represented by a blue star). 
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This intriguing and exciting aspect will be developed elsewhere. 

VI. CONCLUSION. 

This paper gives an example of network dynamics where the nonlinearity induces par- 
ticularly prominent effects that cannot be anticipated by the mere analysis of the graph 
topology. In particular we exhibit a dynamic differentiation in the capacity that a unit has 
to transmit information. We also argue on theoretical grounds, and numerically show (see 
Fig. |2J) that the dynamics differentiation is not revealed by correlation functions. It is purely 
an effect of the dynamics transverse to the chaotic attractor that must be handled with the 
proper tool. We show that the linear response gives quite a bit more information than the 
correlation function, provided that its computation takes into account the singularity of the 
SRB measure transversally to the attractor. This is the case with Ruelle linear response 
theory and this opens the perspective for developing an extension of statistical mechanics 
for the analysis of networks dynamics with nonlinear units. 
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bounded away from zero. Nevertheless we will adopt the point of view defended in |l|. If 
there is a linear response theory for our system, it is necessarily of the form eqs ©-©, since 
there are no reasonable alternative. What could happen is that the sum diverges leading to 
an infinite response. On numerical grounds, one has to check that the time average used to 
compute the ergodic average [see eq. (|12|) ] does not increase with the sample length. 

[13] Note that a linear response theory has also been proposed in 8]. However, it requires the 
invariant measure to have a density. This is only true for the conditional measure along 
unstable manifolds. As a matter of fact, this theory does not contain the stable term. 
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